Mapping CPOL to a cartesian grid

Scott Collis, Argonne National Lab </center>

This notebook is a simple notebook to show the loading, plotting on a map and gridding Darwin CPOL data to a Cartesian grid at two different resolutions


In [27]:
#load relevant modules, set plotting to inline
import pyart
from matplotlib import pyplot as plt
import numpy as np
import wradlib
from pyart.core.echo_classes import *
from matplotlib import rc
from netCDF4 import num2date
print pyart.__version__
%matplotlib inline


1.0.0.dev-293d051

In [48]:
#Data is at: https://engineering.arm.gov/~collis/cpol/

In [3]:
filename = '/data/cpol/Gunn_Pt_rapic_20120315_102002_PPI.UF'
radar = pyart.io.read(filename)

In [5]:
print radar.fields.keys()


['corrected_reflectivity', 'differential_phase', 'cross_correlation_ratio', 'spectrum_width', 'differential_reflectivity', 'velocity']

In [20]:
myd = pyart.graph.RadarMapDisplay(radar)
fig = plt.figure(figsize = [18,10])
myd.plot_ppi_map('corrected_reflectivity', vmin = -16, vmax = 64, resolution = 'h')

m= myd.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary(fill_color='blue')
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 11),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000, 
      runway1y[0]-apannoty+1000, head_width=3000, 
      head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold'})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold'})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold'})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold'})
m.drawmapscale(130.8, -13.2, 131, -13.2, 100, barstyle='fancy', fontcolor='w', fillcolor1='w', fillcolor2='k')


Out[20]:
[<matplotlib.lines.Line2D at 0x118372d10>,
 <matplotlib.lines.Line2D at 0x118384410>,
 <matplotlib.lines.Line2D at 0x118384a50>,
 <matplotlib.lines.Line2D at 0x11838e0d0>,
 <matplotlib.patches.Polygon at 0x11838e710>,
 <matplotlib.patches.Polygon at 0x11838ecd0>,
 <matplotlib.patches.Polygon at 0x11839a1d0>,
 <matplotlib.patches.Polygon at 0x11839a690>,
 <matplotlib.lines.Line2D at 0x11839ab50>,
 <matplotlib.lines.Line2D at 0x118d5f0d0>,
 <matplotlib.lines.Line2D at 0x118d5f710>,
 <matplotlib.text.Text at 0x118d5fd50>,
 <matplotlib.text.Text at 0x118d6a5d0>,
 <matplotlib.text.Text at 0x118d6aa10>,
 <matplotlib.text.Text at 0x118d6ae50>]

In [21]:
cpol_map = pyart.map.grid_from_radars((radar,),  
                                        grid_shape=(35, 401, 401),
                                        grid_limits=((0, 17000), (-150000, 150000), (-150000, 150000)),
                                        grid_origin = (-12.249027,131.044402),
                                        fields=radar.fields.keys(),
                                        refl_field='corrected_reflectivity',
                                        max_refl=100.)

In [22]:
pyart.io.write_grid('/data/cpol/cpol_20120315_102002.nc', cpol_map)

In [47]:
display = pyart.graph.GridMapDisplay(cpol_map)
# create the figure
font = {'size': 14}
rc('font', **font)
fig = plt.figure(figsize=[15, 8])

# panel sizes
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]

# parameters
level = 5
vmin = -16
vmax = 64
lat = -11.89
lon = 130.7

# panel 1, basemap, radar reflectivity a
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution = 'h')
display.plot_grid('corrected_reflectivity', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)

m= display.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 6),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000, 
      runway1y[0]-apannoty+1000, head_width=3000, 
      head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold','size': 12})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold','size': 12})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold', 'size': 12})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold', 'size': 12})
m.drawmapscale(130.8, -13.0, 131, -13.0, 100, barstyle='fancy', fontcolor='b', fillcolor1='b', fillcolor2='k')



# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(cax=cbax)

# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from SGP CF (km)')

# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)

# add a title
slc_height = cpol_map.axes['z_disp']['data'][level]
dts = num2date(cpol_map.axes['time']['data'], cpol_map.axes['time']['units'])
datestr = dts[0].strftime('%H:%M Z on %Y-%m-%d')
title = 'Sliced at ' + str(slc_height) + ' meters at ' + datestr
fig.text(0.5, 0.9, title)


Out[47]:
<matplotlib.text.Text at 0x135d4e490>

In [41]:
cpol_map_hr = pyart.map.grid_from_radars((radar,),  
                                        grid_shape=(35, 601, 601),
                                        grid_limits=((0, 17000), (-150000, 150000), (-150000, 150000)),
                                        grid_origin = (-12.249027,131.044402),
                                        fields=radar.fields.keys(),
                                        refl_field='corrected_reflectivity',
                                        max_refl=100.)

In [44]:
pyart.io.write_grid('/data/cpol/cpol_hr_20120315_102002.nc', cpol_map_hr)

In [46]:
display = pyart.graph.GridMapDisplay(cpol_map_hr)
# create the figure
font = {'size': 14}
rc('font', **font)
fig = plt.figure(figsize=[15, 8])

# panel sizes
map_panel_axes = [0.05, 0.05, .4, .80]
x_cut_panel_axes = [0.55, 0.10, .4, .30]
y_cut_panel_axes = [0.55, 0.50, .4, .30]
colorbar_panel_axes = [0.05, 0.90, .4, .03]

# parameters
level = 5
vmin = -16
vmax = 64
lat = -11.89
lon = 130.7

# panel 1, basemap, radar reflectivity a
ax1 = fig.add_axes(map_panel_axes)
display.plot_basemap(resolution = 'h')
display.plot_grid('corrected_reflectivity', level=level, vmin=vmin, vmax=vmax)
display.plot_crosshairs(lon=lon, lat=lat)

m= display.basemap
runway1x,runway1y=m(np.array([130.865321, 130.89427]),np.array([-12.40945, -12.419184]))
runway2x,runway2y=m(np.array([130.870643, 130.870696]), np.array([-12.409041, -12.422505]))
apannotx, apannoty=m(130.75, -12.37)
armx, army=m(130.891848,-12.424575)
cpolx, cpoly=m(131.044402,-12.249027)
bayx, bayy=m(130.9,-12.3)
m.drawcoastlines()
m.drawmapboundary()
m.drawparallels(np.linspace(-11, -13, 11),labels=[1,0,0,0])
m.drawmeridians(np.linspace(130,132, 6),labels=[0,0,0,1])
m.drawrivers()
m.plot(runway1x, runway1y, 'k-')
m.plot(runway2x, runway2y, 'k-')
# P.arrow( x, y, dx, dy, **kwargs )
#P.arrow( 0.5, 0.8, 0.0, -0.2, fc="k", ec="k",
plt.arrow(apannotx, apannoty, runway1x[0]-apannotx-4000, 
      runway1y[0]-apannoty+1000, head_width=3000, 
      head_length=3000, fc='w', ec='w' )
plt.text(apannotx-5000,apannoty,'Airport', fontdict={'color':'white', 'weight':'bold','size': 12})
m.plot(armx,army,'bo')
plt.text(armx+1500, army-3000, 'ARM Darwin\n site', fontdict={'color':'black', 'weight':'bold','size': 12})
m.plot(cpolx,cpoly,'bo')
plt.text(cpolx+1500, cpoly+1000, 'CPOL Radar', fontdict={'color':'black', 'weight':'bold', 'size': 12})
plt.text(bayx,bayy,'Shoal \n Bay', fontdict={'color':'white', 'weight':'bold', 'size': 12})
m.drawmapscale(130.8, -13.0, 131, -13.0, 100, barstyle='fancy', fontcolor='b', fillcolor1='b', fillcolor2='k')



# colorbar
cbax = fig.add_axes(colorbar_panel_axes)
display.plot_colorbar(cax=cbax)

# panel 2, longitude slice.
ax2 = fig.add_axes(x_cut_panel_axes)
display.plot_longitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)
ax2.set_xlabel('Distance from SGP CF (km)')

# panel 3, latitude slice
ax3 = fig.add_axes(y_cut_panel_axes)
display.plot_latitude_slice('corrected_reflectivity', lon=lon, lat=lat, vmin=vmin, vmax=vmax)

# add a title
slc_height = cpol_map_hr.axes['z_disp']['data'][level]
dts = num2date(cpol_map_hr.axes['time']['data'], cpol_map_hr.axes['time']['units'])
datestr = dts[0].strftime('%H:%M Z on %Y-%m-%d')
title = 'Sliced at ' + str(slc_height) + ' meters at ' + datestr
fig.text(0.5, 0.9, title)


Out[46]:
<matplotlib.text.Text at 0x1349becd0>

In [ ]: